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In the present work, we improve a numerical method, developed to solve the Gross-Pitaevkii non- 
linear Schrodinger equation. A particular scaling is used in the equation, which permits to evaluate 
the wave-function normalization after the numerical solution. We have a two point boundary value 
problem, where the second point is taken at infinity. The differential equation is solved using the 
shooting method and Runge-Kutta integration method, requiring that the asymptotic constants, for 
the function and its derivative, are equal for large distances. In order to obtain fast convergence, 
the secant method is used. 
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Precise and fast numerical solutions to non-linear equa- 
tions have become considerable important in computa- 
tional physics. So, the numerical procedure are rele- 
vant also to be described, when treating such problems, 
considering the computational time consuming. In the 
present work, we pay attention specially to this problem, 
proposing an alternative approach to a method recently 
described in Ref. 0] , which was used to solve the Gross- 
Pitaevkii 1^ nonlinear Schrodinger equation (NLSE) for 
trapped neutral atoms, with positive two-body scatter- 
ing lengths. In Ref. |§, the NLSE treated in Ref. |] was 
extended to a time-dependent one, for both positive and 
negative two-body scattering lengths, where the Crank- 
Nicolson algorithm (appropriate for time evolution) was 
considered. This approach, however, has the disadvan- 
tage that can only reach the stable solutions. In case one 
needs to add other non-linear terms (of higher order) in 
the original equation Q], it is not feasible to reach an- 
other region of stable solutions if in between there is an 
unstable region. This implies that it should be appro- 
priate to combine a static method (as the one used in 
Ref. 1^) with the method used in Ref. |^ when we are 
interested in obtain all the stable and unstable solutions 
and also the corresponding time evolutions. So, in this 
perspective, any improvement of the method considered 
in Ref. ^ would be relevant. 

In the following, we briefly describe the physics re- 
lated to the NLSE considered in 0], and the numerical 
procedure used to solve it. Then we present an alterna- 
tive approach, which can considerably reduce the time to 
search for the solutions and the normalizations. 

The nonlinear Schrodinger equation, which describes 
the condensed wave-function in the mean-field approxi- 
mation can be written as j|] 
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where m is the mass of a single atom, uj the angular 
frequency of the trap, /i the chemical potential and a 
the scattering length. In the present approach, as we are 
more concerned with numerical aspects, for convenience 
we treat only cases with negative scattering lengths (a < 
0) 0. Later we also consider the inclusion of a three body 
interaction term. 

The chemical potential /x is fixed by the number N 
of atoms in the condensed state, which is given by the 
normalization condition: 



(2) 



In Refs. and ||] the NLSE for Bose-Einstein conden- 
sates, as given in Eq. (P, was solved numerically. In Ref. 

, the shooting and the Runge-Kutta methods were 
combined. For a given normalization parameter its was 
solved the corresponding dimensionless equation. The 
asymptotic form of the wave-function was renormalized 
to be equal to the numerical wave-function for a suf- 
ficiently large distance. The wave-function normaliza- 
tion parameter was increased until the Wronskian of the 
asymptotic behavior of the numerical and the analytic 
function changes sign. 

Next, we present in detail the numerical approach we 
have used, in order to show up the similarities and subtle 
differences between this approach and the one of Ref. § 
As we suggest from our experience, such subtle differ- 
ences in the numerical procedures will reduce consider- 
ably the time to search for the solutions. We first rewrite 



*For the numerical considerations, there is no restrictions 
about the sign of the scattering length, as the solutions with 
a > are equally accessible using the same procedure and 
changing the sign of the non- linear term in Eq. (hi). 
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Eq.(|l]) in dimensionless units, in order to become appar- 
ent the physical scales contained in Eq. (|l|) . By rescaling 
Eq.(|l]) for the s-wave solution, we obtain [Q 
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where x = yj2muj/h r, <i>(x) = ^/S^T\a\N r^'(r) and 
/3 = ^/fiuu < 3/2. The normalization for <I>(x), obtained 
from Eqj2[ defines a real number n (given as |C^f | in 
Ref. pi) related to the number of atoms N: 



dx\^{x)\ 



= n, 



where n = A^(87r|a|) 
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We would like to emphasize that, by using this scaling 
procedure, the numerical solutions for the equation are 
free of any normalization constraint, or other parameter 
dependence. The parameter N, related to the number of 
particles, was removed from the differential equation and 
its not necessary to check Eq. (^) or at all steps of 
the calculation. The normalization is found a posteriori, 
using Eq. (Q). 

The boundary conditions for Eq. (^) are given as 



'I>(0)=0 and <^{x)\ 
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where C is a constant to be determined. By using Runge- 
Kutta method and starting with the given $(0), the prob- 
lem is reduced to determine the value of the correspond- 
ing derivative, $'(0), which satisfies the asymptotic con- 
dition at infinity. So, we are tempted to shoot ||^,^ many 
values for $'(0) until we obtain numerically a constant for 
large distances. At a certain Xmax we expect a constant, 
given by 



^numix) exp 
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This process is very laborious and difficult to reach pre- 
cise solutions due to the problem of verifying, for some 
large x, when C$ is constant, within the required numer- 
ical precision. The way to overcome these difficulties is 
to consider the asymptotic derivative of ^{x), that is 
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and also determine (numerically) the expression 
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When we are using the correct value of $'(0) we also 
should obtain C$ = C$' = C for a large enough x — 

Xmax • 

An useful remark we can make is that when we overes- 
timate the value of <i>'(0), increases and C$' decreases. 
The inverse happens when we subestimate $'(0). So, this 
condition is valuable to tune $'(0). It corresponds to 
solve the equation 



(7$ — (7$' = 0, 
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having^' (0) as the unknown variable. Substituting (|^) 
and (g) in (^) we recover the expression for the Wron- 
skian W^($„„m (xmax), $asym (aJmax)) = Stated in [0. 
Eq. (||) can be solved by secant method So, we begin 
with an approximate solution for <&'(0), as an input to 
the secant method, where C$ and C$/ are evaluated by 
the Runge-Kutta method. We should emphasize that, 
to succeed with such method the original guess for $'(0) 
should be not far from the correct value, otherwise the 
method can conduct to the trivial solution ^{x) =0 or 
to overflows. In our procedure, for a fixed /3, Xmax was 
first estimated to be equal to 4.2 and <i>'(0) was used as 
an initial trial to extend Xmax to 5.6 and subsequently to 
7.0. Once we find a solution for $'(0), for a given f3, we 
decrease f3 slightly by A/3 using the previous $'(0) to find 
the new $'(0). This process allows us to "walk" along (3 
values, obtaining the corresponding solutions and results 
for n. 

Although the secant method can become unstable un- 
der certain conditions, in this case it will not occur, as 
we explain in the following. We found that the secant 
method is appropriate, as we can be as near as desired 
to the solution, starting with a given analytical solution 
of the corresponding linear Schrodinger equation. So, we 
just need to implement an automatic algorithm routine 
to make slow variations of /3 and the corresponding slow 
shift (from the initially zero) of ^'(O), in order to satisfy 
the corresponding non-linear equation. In this way, we 
are always near the solution, such that the secant method 
can be applied. We think the same procedure can be gen- 
erally applied for solitonic equations. The algorithm of 
slow variation of /? (deformation algorithm) does not need 
the estimation for the derivative of the wave-function at 
X = 0, as given in Eq.(3.7) of Ref. ||l|) for every solu- 
tion, except for the first one were we take it near the 
harmonic oscillator solution. We understand that the an- 
alytical approximation given in Ref. |l[] to estimate the 
derivative at a; = is not the most convenient in the 
present case. Considering that in general such equations 
are highly non-linear, the initial guess for the derivative 
of $ at x = can easily cause overflow when determining 
the asymptotic constants (Wronskian) at large distances. 
Our initial guess can be very close to the harmonic oscil- 
lator solution, which corresponds to $'(0) = 0, avoiding 
possible overflows for sufficient large distances. In our 
numerical approach, considering that /? = 1.5 is a trivial 
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solution of the linear harmonic oscillator, we started with 
/? = 1.4, trial $'(0) = 0.6 and A/3 = 0.02. For each /3 it 
was necessary 4 ~ 6 iterations in the secant method 
for each Xmax- 

Our results, for several values of /3, are partially listed 
in Table I. The solutions with (3 < 0.4 are unstable 
and not shown in Ref. |3|. However, the solutions with 
/3 > 0.4 well agree with their results. As one can observe 
in Table I, the method also can reach solutions with neg- 
ative chemical potentials (/3 < 0). A numerical stability 
check, which can be done by evolving the static solutions, 
can easily be followed by using a time dependent method, 
as the Crank-Nicolson method . 

In Fig. 1 we also show three plots for the chemical 
potential, (3, as a function of rt, in case of zero angular 
momentum. The three plots shown correspond to the 
lower radial states (n^ = 0,2,4) of Eq. (|^). The plot 
labeled with = corresponds to Table L In the limit 
of the harmonic oscillator solution, where n = and the 
equation is linear, we obtain the usual known solutions. 

In Fig. 2 we have another example of the application 
of the method described here. In this case, we consider 
the addition of another non-linear term inside the square 
brackets of Eq. (||) , given by 
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which can be directly related to the three-body effects, 
where gs is the non-dimensional strength of the corre- 
sponding three-body interaction. The physical conse- 
quences of the addition of such a term in the NLSE is 
discussed in both references given in 



TABLE I. Numerical solutions 
two-body interaction, when the two- 
ative. We consider that at Smax ~ 
asymptotic limit. 



for the NLSE including 
-body scattering a is neg- 
7.0 we have achieved the 
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FIG. 1. We show the chemical potential, /3, as a function 
of n, which is related to the number of particles by the Eq.(^. 
The three plots shown correspond to the lower radial states 
(ur) of Eq. (^) (with zero angular momentum), such that in 
the limit of the harmonic oscillator solution, where n = 0, we 
have the usual known solutions. The plot labeled with = 
corresponds to Table I. 
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FIG. 2. We show the chemical potential, /3, as a function 
of n, in case we consider three-body effects. The three plots 
shown correspond to the radial state solutions for gs = Q (no 
quintic term), gz = 0.016 and gz = 0.03. The plot labeled 
with Qs = G corresponds to Table I. 

To finalize, we have presented in detail an improvement 
to a numerical procedure used to solve a non-linear dif- 
ferential equation, which is commonly applied for Bose- 
condensed states. In our example, we solve the Gross- 
Pitaevskii equation with an attractive two-body and a 
repulsive three-body interaction. We should note that, 
by using a simplified scaling procedure [given in Eqs.(||) 
and (^J, the numerical solutions for the equation are free 
of any normalization constraint, or other parameter de- 
pendence. The parameter iV, related to the number of 
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particles, was removed from the differential equation and 
it is not necessary to obtain the normalization at each 
step of the calculation. Eq. gives the normalization 
a posteriori. So, by using the above scaling procedure, it 
emegers the main difference between the present method 
and the one given in Ref. 0, when looking for the solu- 
tions of the NLSE: 

- in our approach, we searched for the derivative of 
the wave-function at a; = till the asymptotic con- 
stants match (when the Wronskian vanishes), and 
the normalization is given at the end; 

- in 01 the normalization parameter A is incremented 
till the sign of the Wronskian is reached. For the 
final renormalization they also use other interme- 
diate parameters, as Aq, iVo, Ai, Ni. 

As we are not restricted by the normalization, our ap- 
proach is a clear improvement to the method given in 
Ref. [Q , particularly when considering the simplification 
and the transparency in the normalization procedure. 
Such advantage can be further explored when more in- 
volved calculation are presented, as in Ref. [7] where 
collective excitations are evaluated. A different scaling 
removes the normalization constraint and allows one to 
obtain it a posteriori, after the numerical solutions are 
achieved for the eigenvalues. 

In our numerical procedure, we employ the shoot- 
ing method on Runge-Kutta integration, matching the 
asymptotic constants for the wave-function and for the 
corresponding derivative. This procedure is shown to be 
equivalent to make the Wronskian vanish at such large 
distances. In order to obtain a faster convergence to 
the solution, we also included the secant method. The 
numerical optimization to the method employed in Ref. 
0, described here, is not restricted to the NLSE we have 
used. It can be used quite generically for second order 
solitonic differential equations whose solutions asymptot- 
ically vanish at large distances. We consider particularly 
relevant an optimization of the method of Ref. [Q in the 
perspective of looking for solutions of differential equa- 
tions with higher order non-linear terms; and also in the 
perspective of combining such method (appropriate for 
static solutions) with a time-dependent one. 
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